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LONG-TERM  GOALS 

To  create  a  suite  of  computer  codes  (programs)  that  accurately  model  wave  propagation  and  scattering 
in  shallow  water,  to  quantify  the  effects  of  poroelastic,  stratified  sediment  and  elastic  mode  conversion 
on  wave  propagation  in  shallow  water  at  frequencies  ranging  from  10  Hz  to  several  hundred  kHz. 

To  relate  these  effects  to  sonar  system  behavior  and  efficient  design  of  sonar  capable  of  finding 
buried/semi-buried/proud  mines  in  shallow  water  environments,  in  the  presence  of  poroelastic 
sediment  and  bubble  layers. 

To  understand  the  fundamental  linear  and  nonlinear  physics  of  wave  propagation  in  poroelastic 
sediment  with  bubbles. 

OBJECTIVES 

To  continue  development  and  testing  of  K-space  Integral  Equation  codes  (K-SIE),  finite  difference 
time  domain  codes  (FDTD),  and  Hybrid  Parabolic  marching/Finite  difference  (HP-MFD)  and  Hybrid 
Integral  Equation/Finite  Difference  (HIEFD)  methods  for  fast  and  accurate  scattering  and  inverse 
scattering  of  objects  buried  in  sediment  with  poroelastic  (Biot)  properties  and  containing  bubble 
populations. 

To  develop  a  Java-based,  interactive,  web  site  http://borup.aim.utah.edu  [41]  to  make  pseudo-analytic, 
FDTD,  Integral  Equation,  and  HIEFD  and  HP-MFD  calculation  methods  in  3D  and  2D,  available  to 
other  researchers  over  the  web. 

RESULTS 

Our  results  encompass  five  main  types  of  codes:  (1)  FDTD  (Finite  Difference  Time  Domain);  (2)  K- 
SIE  (K-Space  Integral  Equation);  (3)  P-AC  (pseudo-analytic  code);  and  the  two  hybrid  methods:  (4) 


Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 
VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 
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HIEFD  (Hybrid  Integral  Equation  Finite  Difference);  and  (5)  HP-MFD  (Hybrid  Parabolic-Marching 
Finite  Difference).  Concepts  found  in  references  [1-45]  have  proven  to  be  very  useful  in  this  research. 

1.  Results  with  FDTD  codes  (Finite  Difference  Time  Domain) 

We  have  finished  the  initial  phases  of  the  adaptation  of  a  FDTD.  2D,  geophysics  Biot  code  (given  to 
us)  to  the  P-SV  wave  elastic  Biot  [12]  problem  in  2D.  In  doing  so  we  uncovered  an  error  in  the  gifted 
code,  and  have  generalized  it  to  include  non-homogeneities  in  porosity,  permeability,  bulk  moduli, 
shear  modulus,  and  viscosity  of  the  pore-filling  fluid. 

We  have  incorporated  our  “pill-box”  scheme  into  a  hybrid  IEFD  for  the  acoustic  scalar  case,  wherein, 
the  IE  code  is  used  outside  and  the  FDTD  code  is  used  inside  the  same  pillbox. 

Forward  scattering  code  for  the  time  domain  scattering  in  two  and  three  dimensions  with  the  advanced 
Berenger-type  Perfectly  Matched  Fayer  absorbing  boundary.  These  effective  boundary  conditions 
reduce  the  anomalous  scattering  from  grid  edges  by  71  dB,  versus  approximately  40  dB  reduction 
provided  by  the  standard  absorbing  BC's  of  Mur  [31]. 

We  have  extended  the  2D  and  3D  forward  scattering  time  domain  code  to  include  the  effects  of 
dispersion  as  dictated  by  the  Kramers- Kronig  relations  in  two  distinct  ways:  (1)  using  a  Z-transform 
approach  developed  by  Sullivan  [43];  and  (2)  a  state  variable  approach  developed  by  us  [31].  These 
alternative  methods  are  compared  to  one  another  for  accuracy.  Furthermore,  we  have  incorporated  the 
new  PMF  absorbing  BC's  in  the  [31]  state  variable  approach. 

2.  Results  with  K-SIE  codes  (K-Space  Integral  Equation) 

a.  Dictionary.  By  encoded  the  Biot  Green’s  function  (point  response  function)  as  described  in  [Boutin 
et  al  we  have  uncovered  an  important  physical  misinterpretation  that  could  render  the  point  response 
function,  as  displayed  in  [8]  incomplete  from  a  physical  point  of  view.  This  is  true  even  though  on  a 
purely  mathematical  level  the  Green’s  function  is  strictly  correct.  Furthermore,  we  have  established  a 
convenient  dictionary  that  allows  one  to  pass  from  the  set  of  terminology  favored  by  the  “French 
school”  (Boutin,  Bourbie’  Bonnet,  etc.)  [8,18,19,28]  to  the  terminology  common  to  Stoll  [29], 
Kibblewhite  [30],  Johnson  and  Plona  [1],  etc.  This  discrepancy  in  terminology  arises  primarily  from 
the  use  of  homogenization  theory  in  the  French  school,  and  a  stricter  adherence  to  the  Biot  approach  in 
Stoll/Kibblewhite  [29,30]. 

We  are  also  in  the  midst  of  carrying  out  model  wave  propagation  simulations  with  the  parameters 
espoused  by  Chotiros  [39].  We  have  generalized  the  method  of  Kibblewhite  &  Wu  [30]  for  the 
determination  of  reflection  coefficients  in  visco-elastic  porous  layered  half-spaces  to  yield  a  layered 
Biot  point  response  function  for  arbitrary  layered  porous  media,  using  the  technique  we  previous 
developed  for  the  acoustic  case  [42]. 

We  have  developed  a  dictionary  for  the  various  notations  employed  by: 

“The  Rock  Physics  Handbook”,  Mavko  et  al.  [27], 

“Wave  Interactions  as  a  Seismo-acoustic  Source”,  Kibblewhite  and  Wu  [30] 

“Acoustics  of  Porous  Media”,  Bourbie’  et.  al.  [28] 

“Sediment  Acoustics”,  Robert  Stoll  [29], 

Several  papers  (Johnson  &  Plona  [1],  Boutin  et.  al.[8],  Bonnet  et.  al.  [8,18,19],  Chotiros)  [39]. 


An  example  of  the  relationship  and  utility  of  such  cross-reference  is  the  following:  Defining 
P=porosity,  u=averaged  “solid”  displacement,  U=averaged  “fluid”  displacement,  w=|3(u-U)  the 
“increment  of  fluid  content”  of  [Stoll,  1989],  we  can  write: 

//V2u  +(//  -  2 p)V e  -  CVg  =  pun  -  Pf\\tt 
CVe-  MVc  =  p f  u tt  -  mwtt  -  (///w)w; 

In  the  above  equations,  e,  p,  H,  C,  M,  p,pf,  rj,  k  are  all  as  defined  by  Stoll  [29].  These  equations  can 
be  rewritten  in  the  form  used  the  "French  school"  [8,18,19].  This  results  in  the  identification  of 
parameters:  for  example,  the  ratio  Cl M  =  aBoutin  =  1-  Kb  /  Ks  in  the  Boutin  [8],  which  in  turn  is 

essentially  the  "Biot  coefficient"  discussed  in  Bourbie's  book  [28].  Also  |3  serves  as  porosity  in 
Stoll/Kibblewhite  notation,  but  is  (3  =  (a  -  <ft)l  Ks  +  <j) / Kj  in  the  Boutin  et  al.  [8],  although  Boutin  uses 

n  not  (j>  for  porosity.  This  |3  is  equivalent  to  1/M  in  the  notation  of  Bourbie'  [28] 

A  careful  analysis  of  the  equations  reveals  that  the  m  parameter  in  Stoll/Kibblewhite  notation  m  = 
ipf/(|),  can  be  replaced  by  m  — »  m  +  qF2/(Kco).  where  F2  is  the  imaginary  part  of  the  complex  F,  the 
Biot  frequency  dependent  structure  factor,  and  plays  a  role  analogous  to  the  imaginary  part  of  H,  the 
inverse  of  the  generalized  Darcy  coefficient  of  homogenization  theory.  The  real  part  of  F,  FR  is 
analogous  to  the  real  part  of  the  inverse  H  of  the  generalized  Darcy  coefficient  of  homogenization 
theory. 

These  identifications  allow  for  accurate  model  parameters  (with  physically  relevant  values)  to  be  input 
into  the  Biot  Green's  function  developed  within  the  context  of  homogenization  theory.  These 
identifications  also  allow  for  theory  to  be  compared  with  experiment  and  thus  ultimately  inversion 
methods  based  on  reduction  of  residual  (model  minus  measurement)  error. 

b.  Free  space  background  codes.  The  2D  fom’ard  and  inverse  code  for  the  scattering  of  an  elastic 
object  in  an  elastic  background  medium,  with  Lame'  parameter  X  allowed  to  vary,  and  p  held  fixed  is 
running.  The  2D  forward  scattering  code  is  running  for  the  scattering  of  an  elastic  wave  off  of  an 
elastic  object  in  a  fully  elastic  medium:  i.e.  both  X  and  p  are  allowed  to  vary.  The  3D  forward 
scattering  code ,  for  elastic  scattering  from  an  elastic  inhomogeneity  within  an  elastic  medium  (Lame' 
parameter  X  is  spatially  dependent,  but  p  is  held  to  be  fixed  spatially)  has  been  written  and  is  being 
tested  against  certain  pseudo-analytic  solutions.  We  have  also  finished  the  development  and  testing  of 
the  associated  3D  inversion  code  in  transmission  mode,  and  are  in  the  process  of  generalizing  to  the 
important  problem  of  multi-static  reflection  mode  inversion.  We  have  improved  the  convergence  of 
our  K-SIE  codes  by  replacing  the  forward  problem  solver  with  the  Stabilized  version  of  Bi-CG  [34]. 
We  have  also  improved  the  overall  convergence  of  the  inverse  problem  by  use  of  the  fully  nonlinear 
Ribiere-Polak  conjugate  gradient  method. 

The  elastic  Green's  function  in  the  free  space  case  has  a  well-known  non-integrable  singularity  which 
we  handle  by  careful  use  of  sine  functions,  sinc(x)  =  [sin(7tx)]/[7tx],  and  sinc-convolution.  We  also 
address  the  problem  of  the  existence  of  derivative  operators  acting  inside  the  integral  on  the 
displacement  functions  by  means  of  augmenting  the  integral  equations  of  motion  with  explicit  integral 
equations  for  the  strain  tensor  components. 

c.  Born  analysis  of  the  elastic  imaging  problem  We  have  also  carried  out  a  theoretical  analysis  based 
upon  the  Born  assumption,  and  the  integral  equation  formulation  for  elastic  scattering,  that  determines 
what  set  of  sediment  elastic  parameters  (X,  p,  p)  can  be  successfully  inverted  when  sources  and 


receivers  are  located  in  water.  We  derived  the  Born  approximation  to  the  volume  elastic  integral 
equation  for  an  elastic  object  imbedded  in  an  elastic  medium  in  order  to  examine  the  invertability  for 
the  three  isotropic  media  parameters,  p,  p,  X.  The  analysis  shows  that  all  three  parameters  can  be 
imaged  independently  if  a  shear  wave  source  is  available.  Also,  only  a  pressure  receiver  is  needed  (s- 
source,  s-receiver  data  is  not  needed).  In  2D,  it  turns  out  that  p  can  be  independently  imaged  with  a 
shear  source  and  pressure  receiver  (i.e.,  only  perturbations  in  p  scatter  for  this  case  in  2D).  The  other 
two  parameters,  p  and  X,  can  then  be  imaged  with  p-source,  p-receiver  data  in  an  analogous  manner  to 
p  and  k  for  the  acoustic  imaging  problem. 

d.  Layered  Background  Code.  We  have  developed  a  layered  Green's  function  and  concomitant  K- 
space  integral  equation  for  arbitrarily  layered  media.  The  construction  we  have  employed  mirrors  the 
construction  in  [42],  and  is  specifically  designed  to  preserve  the  use  of  the  Fast  Fourier  Transform 
(FFT)  even  in  the  vertical  direction  (the  direction  of  the  arbitrary  layering). 

These  2D  and  3D  codes  and  Green's  functions  are  used  in  two  contexts:  (1)  volume  scattering  for 
remote  sources  and  receivers,  and  (2)  as  a  propagator  within  the  context  of  the  hybrid  FDIE  (finite 
difference/integral  equation)  method.  When  the  codes  are  utilized  in  the  context  of  volume  scattering 
special  care  must  be  taken  to  handle  the  singular  nature  of  the  self-term.  We  have  utilized  "sine" 
functions,  and  a  unique  expression  for  the  sine-convolved  Green's  function  (a  pseudo-analytic  series 
expansion  in  orthogonal  polynomials),  which  was  developed  at  CIPIT,  to  address  this  problem. 

3.  Results  with  P-AC  (pseudo-analytic  code) 

We  have  developed  so-call  Pseudo-analytic  solutions.  These  solutions,  based  on  Bessel  and  Hankel 
function  expansions,  are  embodied  in  code  developed  at  our  Center  for  Inverse  Problems,  Imaging  and 
Tomography  (CIPIT),  and  will  be  made  available  over  the  Web,  in  a  manner  similar  to  the  interface 
that  is  presently  available  over  the  net  at  http://borup.aim.utah.edu  .  This  code  will  enable  researchers 
from  remote  sites  to  test  their  developed  code  (finite  difference,  marching  methods,  boundary  integral 
equations,  etc.)  against  analytic  solutions  for  simple  objects  (layered  spheres  and  cylinders,  and  finite 
groups  thereof)  in  free  space.  For  a  single,  layered,  analytic  object,  there  is  no  inversion,  only 
calculation.  For  two  such  objects,  there  is  only  a  recursion  (even  in  2D  or  3D).  Extension  to  more  or 
overlapping  objects  is  possible. 

4.  Results  with  HIEFD  codes  (Hybrid  Integral  Equation  Finite  Difference) 

The  2D  pillbox  code  has  been  written  and  tested;  further  testing  is  underway.  See  Figure  1  for  a 
geometrical  description  of  the  location  of  the  pillbox  with  respect  to  the  ocean  and  sediment. 

The  2D  "pill-box"  approach  requires  that  equivalent  sources  be  located  on  the  boundary  of  a  finite 
difference  grid  to  allow  the  external  incident  field  to  be  propagated  into  the  grid,  where  it  can  interact 
with  field  inhomogeneities  and  rescatter.  The  external  field  is  propagated  to  the  pillbox  in  this  case  by 
the  fast  calculation  of  the  layered  impulse  response  function  (a  Hankel  transform  is  evaluated).  Since 
the  layered  impulse  response  function  is  a  frequency  domain  object  the  field  is  calculated  at  several 
frequencies,  then  a  Fast  Fourier  Transform  (FFT)  is  applied  to  transform  the  field  on  the  boundary  to 
the  time  domain  (the  equivalent  sources  mentioned  above).  When  the  interaction  has  been  completely 
simulated  via  finite  differences  (FD),  the  field  is  then  re-transformed  back  to  the  Fourier  domain,  and 
the  Green's  function,  and  its  normal  derivative,  are  used  to  propagate  the  fields  back  to  receivers. 

The  finite  difference  codes  required  to  complete  the  Hybrid  IEFD  method  have  been 
completed  and  checked  in  2D  and  3D 


Figure  1 


Cutaway  view  of  ocean,  sediment  and  buried  pillbox  as 
employed  in  the  hybrid  pillbox  method. 


5.  Results  with  HP-MFD  codes  (Hybrid  Parabolic-Marching  Finite  Difference)The  2D  and  3D  pillbox 
propagators  have  been  coded,  and  tested.  The  2  and  3D  pillbox  propagators  use  parabolic-marching 
methods  for  increased  speed  over  HIEFD  methods,  but  at  the  price  of  some  loss  of  accuracy.  The  fast 
calculation  feature  of  the  parabolic-marching  method  propagates  the  external  field  to  the  pillbox  in  this 
case.  After  matching  boundary  conditions  at  the  pillbox  surface,  the  field  is  then  propagated  by  FDTD 
methods  inside  of  the  pillbox.  The  scattered  field  on  the  surface  of  the  pillbox  is  matched  to  the 
parabolic-marching  method  and  then  propagated  through  the  sediment  to  the  ocean  and  then  outward 
to  required  receivers.  In  the  case  of  a  homogeneous  ocean,  calculation  time  may  be  saved  by  use  of 
Green’s  functions  for  in-ocean  propagation,  for  both  hybrid  codes.  See  Figure  1  for  pillbox  geometry. 

IMPACT/APPLICATIONS 

The  development  of  the  generalized  layered  elastic  porous  Green’s  function  mentioned  above,  will 
enable,  with  the  combination  of  this  approach  with  the  FDTD  code  mentioned  earlier,  for  the  2D 
elastic  porous  media,  the  modeling  of  reverberation  and  scattering  from  buried  and  semi-buried  objects 
in  sediment.  The  hybrid  approach  enables  the  efficient  modeling  of  high  contrast,  as  well  as  low 
contrast  buried  mines,  as  insonified  from  remote  sources  and  receivers.  The  use  of  the  Green’s 
propagator  allows  the  accurate  propagation  from  remote  receivers  into  the  porous  sediment,  which  may 
vary  horizontally  in  porosity,  permeability,  bulk  moduli,  shear  moduli,  viscosity,  etc. 


The  Navy  is  very  interested  in  applying  recently  developed  remotely  operated  vehicle  (ROV’s) 
technology  to  the  problem  of  mine  detection  via  multi-static  sonar.  The  hybrid  methods  developed  by 
our  group  and  made  available  over  the  Web  will  enable  the  accurate  and  efficient  modeling  of  these 
multi- static,  multi-platform  sonar  experiments. 

The  K-SIE  methods  allow  the  accurate  propagation  over  large  distances,  and  the  PE  methods  allow 
propagation  over  range  dependent  environments.  Our  PE  code  will  propagate  to  the  buried  or  semi- 
buried  pillbox  in  the  range  dependent  case,  and  the  K-SIE  methods  will  be  used  in  the  range- 
independent  environment.  Once  at  the  pillbox,  whichever  method  is  used,  the  FDTD  codes  will  be 
employed  to  determine  the  accurate  scattered  field  even  for  high  contrast  objects  that  are  buried  in 
porous  elastic  sediment.  The  scattered  field  is  then  Fourier  transformed  at  the  boundary  of  the  pillbox, 
and  propagated  to  remote  receivers  either  in  mult-static  or  mono-static  configurations. 

For  extremely  high  contrast  objects,  the  method  recently  developed  by  Ghosh  Roy  and  co-workers  will 
enable  the  detection  of  buried  and  semi-buried  mines  in  sediment  (i.e.  high  contrast  with  density,  speed 
of  sound,  etc.).  It  is  particularly  well  suited  to  the  multistatic  case.  For  monostatic  data  from  high 
contrast  objects,  the  FDTD  forward  code  and  inversion  is  most  appropriate.  The  forward  scattering 
time  domain  code  employs  4th  order  differencing  with  Berenger  PMF  [31[  absorbing  boundary 
conditions  that  enable  the  accurate  simulation  of  larger  sized  areas  of  seabed  The  FDTD  forward  code 
exists  in  2D  and  3D  forms.  The  associated  inversion  code  exists  in  the  2D  version,  and  employs  the 
Tarantola  [44,45]  method  of  time  domain  regularized  inversion.  It  has  been  tested,  both  the  time 
domain  calculation  of  the  Jacobian  and  its  adjoint  have  been  coded  and  tested. 

By  extending  the  size,  accuracy  and  efficiency  of  our  codes  to  more  accurately  reflect  recent  advances 
in  understanding  sediments  and  Biot  theory  we  increase  the  capability  of  the  Navy  to  find  buried/semi- 
buried  mines  in  surf- zones.  In  particular,  in  keeping  with  our  bio-engineering  background  we  desire  to 
mimic  the  remarkable,  documented  ability  of  the  dolphin  to  detect  prey  and  other  objects  buried  in 
sediment  using  frequencies  ranging  from  10  to  150  kHz  [39]. 

We  will  also  combine  the  results  of  our  research  in  the  use  of  low  frequency  EM  (diffusion 
approximation)  using  some  variant  of  Dempster-Schaeffer  evidential  theory,  and  standard  techniques 
of  data  fusion  to  greatly  enhance  the  capability  of  environmental  scientists,  and  navy  researchers  to 
detect  and  identify  buried  ordnance  in  bubble-containing  and  reverberant  surf- zones. 

We  are  also  contributing  to  the  fundamental  understanding  of  wave  scattering,  transmission  and 
reverberation  in  shallow  water,  and  at  low,  medium  and  high  frequencies  by  being  able  to  compare 
data  collected  in  the  ocean  environment  against  our  codes,  which  accurately  model,  attenuation, 
Kramers- Kro nig  type  dispersion,  poro-elastic  (Biot)  effects,  and  stratification. 

We  have  also  contributed  to  theoretical  understandings  by  explicitly  relating  the  correspondence 
between  the  various  parameters  used  in  homogenization  theory  approaches  adopted  by  the  "French 
school"  (Bourbie1,  Boutin,  Bonnet,  Bard,  Kazi-Aoul,  etc.)  [8,18,19,28]  and  the  Biot  parameters  as 
defined  originally  by  Biot,  and  adapted  by  Chotiros,  Stoll,  Kibblewhite  and  Wu,  Stern,  and  others 
[39,29,30].  In  particular  the  generalized  complex  Darcy  coefficient  arising  in  Boutin  et  al.  [8],  is 
related  to  the  equivalent  mass  parameter  of  Stoll  and  others,  and  Biot  theory. 


TRANSITIONS 


Transitioning  to  others,  of  technology  developed  under  this  project  and  funding,  are  of  two  types:  (1) 
use  of  this  technology  within  our  group  and  industrial  partners  (e.g.:  TechniScan,  Inc  [36],  EMI  [37]; 
and  (2)  use  of  this  technology  outside  of  this  larger  set.  In  the  first  category,  we  and  our  sponsors  [38] 
for  breast  cancer  imaging  and  detection,  are  benefiting  from  advances  in  parabolic  methods  for  fast 
inverse  imaging  of  speed  of  sound,  which  is  one  prominent  indicator  for  cancer.  Also,  the  AFOSR 
sponsored  project  for  detecting  buried  objects  by  ground  penetrating  radar,  has  benefited  by  use  of  the 
PML  FDTD  codes.  In  the  second  category,  our  PMF  FDTD  codes  have  enhanced  the  work  of  EMI 
[37]  and  Waag  [40], 

RELATED  PROJECTS 

Some  the  related  projects  are  described  in  the  first  category  listed  above  in  the  TRANSITONS  section. 
A  further  related  project  is  sponsored  by  the  ONR  and  used  EM  Diffusion  imaging  to  detect  buried 
mines  in  the  surf  zone  (J.  Kravitz,  Ph.D.,  program  manager).  The  Green’s  function  imaging  theory  can 
be  used  directly  or  in  an  approximate,  but  faster  form  called  EM  Migration. 
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